home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 3: Developer Tools / Linux Cubed Series 3 - Developer Tools.iso / devel / lang / lisp / guile-ii.src / guile-ii / guile-src / slib / prime.scm < prev    next >
Encoding:
Text File  |  1994-05-25  |  5.7 KB  |  150 lines
  1. ;;;; "prime.scm", prime test and factorization for Scheme
  2. ;;; Copyright (C) 1991, 1992, 1993 Aubrey Jaffer.
  3. ;
  4. ;Permission to copy this software, to redistribute it, and to use it
  5. ;for any purpose is granted, subject to the following restrictions and
  6. ;understandings.
  7. ;
  8. ;1.  Any copy made of this software must include this copyright notice
  9. ;in full.
  10. ;
  11. ;2.  I have made no warrantee or representation that the operation of
  12. ;this software will be error-free, and I am under no obligation to
  13. ;provide any services, by way of maintenance, update, or otherwise.
  14. ;
  15. ;3.  In conjunction with products arising from the use of this
  16. ;material, there shall be no use of my name in any advertising,
  17. ;promotional, or sales literature without prior written consent in
  18. ;each case.
  19.  
  20. (require 'random)
  21. (require 'modular)
  22.  
  23. ;;; (modulo p 16) is because we care only about the low order bits.
  24. ;;; The odd? tests are inline of (expt -1 ...)
  25.  
  26. (define (prime:jacobi-symbol p q)
  27.   (cond ((zero? p) 0)
  28.     ((= 1 p) 1)
  29.     ((odd? p)
  30.      (if (odd? (quotient (* (- (modulo p 16) 1) (- q 1)) 4))
  31.          (- (prime:jacobi-symbol (modulo q p) p))
  32.          (prime:jacobi-symbol (modulo q p) p)))
  33.     (else
  34.      (let ((qq (modulo q 16)))
  35.        (if (odd? (quotient (- (* qq qq) 1) 8))
  36.            (- (prime:jacobi-symbol (quotient p 2) q))
  37.            (prime:jacobi-symbol (quotient p 2) q))))))
  38.  
  39. ;;;; Solovay-Strassen Prime Test
  40. ;;;   if n is prime, then J(a,n) is congruent mod n to a**((n-1)/2)
  41.  
  42. ;;; See:
  43. ;;;   Robert Solovay and Volker Strassen,
  44. ;;;     "A Fast Monte-Carlo Test for Primality,"
  45. ;;;     SIAM Journal on Computing, 1977, pp 84-85.
  46.  
  47. ;;; checks if n is prime.  Returns #f if not prime. #t if (probably) prime.
  48. ;;;   probability of a mistake = (expt 2 (- prime:trials))
  49. ;;;     choosing prime:trials=30 should be enough
  50. (define prime:trials 30)
  51. ;;; prime:product is a product of small primes.
  52. (define prime:product
  53.   (let ((p 210))
  54.     (for-each (lambda (s) (set! p (or (string->number s) p)))
  55.       '("2310" "30030" "510510" "9699690" "223092870"
  56.            "6469693230" "200560490130"))
  57.     p))
  58.  
  59. (define (prime:prime? n)
  60.   (set! n (abs n))
  61.   (cond ((<= n 36) (and (memv n '(2 3 5 7 11 13 17 19 23 29 31)) #t))
  62.     ((= 1 (gcd n prime:product))
  63.      (do ((i prime:trials (- i 1))
  64.           (a (+ 1 (random (- n 1))) (+ 1 (random (- n 1)))))
  65.          ((not (and (positive? i)
  66.             (= (gcd a n) 1)
  67.             (= (modulo (prime:jacobi-symbol a n) n)
  68.                (modular:expt n a (quotient (- n 1) 2)))))
  69.           (if (positive? i) #f #t))))
  70.     (else #f)))
  71.  
  72. ;;;;Lankinen's recursive factoring algorithm:
  73. ;From: ld231782@longs.LANCE.ColoState.EDU (L. Detweiler)
  74.  
  75. ;                  |  undefined if n<0,
  76. ;                  |  (u,v) if n=0,
  77. ;Let f(u,v,b,n) := | [otherwise]
  78. ;                  |  f(u+b,v,2b,(n-v)/2) or f(u,v+b,2b,(n-u)/2) if n odd
  79. ;                  |  f(u,v,2b,n/2) or f(u+b,v+b,2b,(n-u-v-b)/2) if n even
  80.  
  81. ;Thm: f(1,1,2,(m-1)/2) = (p,q) iff pq=m for odd m.
  82.  
  83. ;It may be illuminating to consider the relation of the Lankinen function in
  84. ;a `computational hierarchy' of other factoring functions.*  Assumptions are
  85. ;made herein on the basis of conventional digital (binary) computers.  Also,
  86. ;complexity orders are given for the worst case scenarios (when the number to
  87. ;be factored is prime).  However, all algorithms would probably perform to
  88. ;the same constant multiple of the given orders for complete composite
  89. ;factorizations.
  90.  
  91. ;Thm: Eratosthenes' Sieve is very roughtly O(ln(n)/n) in time and
  92. ;     O(n*log2(n)) in space.
  93. ;Pf: It works with all prime factors less than n (about ln(n)/n by the prime
  94. ;    number thm), requiring an array of size proportional to n with log2(n)
  95. ;    space for each entry.
  96.  
  97. ;Thm: `Odd factors' is O((sqrt(n)/2)*log2(n)) in time and O(log2(n)) in
  98. ;     space.
  99. ;Pf: It tests all odd factors less than the square root of n (about
  100. ;    sqrt(n)/2), with log2(n) time for each division.  It requires only
  101. ;    log2(n) space for the number and divisors.
  102.  
  103. ;Thm: Lankinen's algorithm is O(sqrt(n)/2) in time and O((sqrt(n)/2)*log2(n))
  104. ;     in space.
  105. ;Pf: The algorithm is easily modified to seach only for factors p<q for all
  106. ;    pq=m.  Then the recursive call tree forms a geometric progression
  107. ;    starting at one, and doubling until reaching sqrt(n)/2, or a length of
  108. ;    log2(sqrt(n)/2).  From the formula for a geometric progression, there is
  109. ;    a total of about 2^log2(sqrt(n)/2) = sqrt(n)/2 calls.  Assuming that
  110. ;    addition, subtraction, comparison, and multiplication/division by two
  111. ;    occur in constant time, this implies O(sqrt(n)/2) time and a
  112. ;    O((sqrt(n)/2)*log2(n)) requirement of stack space.
  113.  
  114. (define (prime:f u v b n)
  115.   (if (<= n 0)
  116.       (cond ((negative? n) #f)
  117.         ((= u 1) #f)
  118.         ((= v 1) #f)
  119.         ; Do both of these factors need to be factored?
  120.         (else (append (or (prime:f 1 1 2 (quotient (- u 1) 2))
  121.                   (list u))
  122.               (or (prime:f 1 1 2 (quotient (- v 1) 2))
  123.                   (list v)))))
  124.       (if (even? n)
  125.       (or (prime:f u v (+ b b) (quotient n 2))
  126.           (prime:f (+ u b) (+ v b) (+ b b) (quotient (- n (+ u v b)) 2)))
  127.       (or (prime:f (+ u b) v (+ b b) (quotient (- n v) 2))
  128.           (prime:f u (+ v b) (+ b b) (quotient (- n u) 2))))))
  129.  
  130. (define (prime:factor m)
  131.   (if
  132.    (negative? m) (cons -1 (prime:factor (- m)))
  133.    (let* ((s (gcd m prime:product))
  134.       (r (quotient m s)))
  135.      (if (even? s)
  136.      (append
  137.       (if (= 1 r) '() (prime:factor r))
  138.       (cons 2 (let ((s/2 (quotient s 2)))
  139.             (if (= s/2 1) '()
  140.             (or (prime:f 1 1 2 (quotient (- s/2 1) 2))
  141.                 (list s/2))))))
  142.      (if (= 1 s) (or (prime:f 1 1 2 (quotient (- m 1) 2)) (list m))
  143.          (append (if (= 1 r) '()
  144.              (or (prime:f 1 1 2 (quotient (- r 1) 2)) (list r)))
  145.              (or (prime:f 1 1 2 (quotient (- s 1) 2)) (list s))))))))
  146.  
  147. (define jacobi-symbol prime:jacobi-symbol)
  148. (define prime? prime:prime?)
  149. (define factor prime:factor)
  150.